﻿using LatoolNet;

namespace NumericalSolution {
  class Poisson1DFV {

    private int m_numberOfControlVolumes;
    private double m_deltax;
    private double m_beta;

    public int NumberOfControlVolumes {
      get { return m_numberOfControlVolumes; }
      set {
        m_numberOfControlVolumes = value;
        m_deltax = 1.0 / m_numberOfControlVolumes;
      }
    }

    public double Beta {
      get { return m_beta; }
      set { m_beta = value; }
    }

    public double DeltaX {
      get { return m_deltax; }
    }

    public Matrix Solve() {

      int bandWidth = 3;
      Matrix mat = new Matrix(m_numberOfControlVolumes, 
                              m_numberOfControlVolumes, bandWidth);
      Matrix vec = new Matrix(m_numberOfControlVolumes, 1);

      double T_L = 0.0;
      double T_R = 1.0;

      mat[0, 0] = 3 / m_deltax;
      mat[0, 1] = -1 / m_deltax;
      vec[0, 0] = m_beta * m_deltax + 2 * T_L / m_deltax;

      for (int i = 1; i < m_numberOfControlVolumes - 1; i++) {
        mat[i, i - 1] = -1 / m_deltax;
        mat[i, i] = 2 / m_deltax;
        mat[i, i + 1] = -1 / m_deltax;
        vec[i, 0] = m_beta * m_deltax;
      }

      mat[m_numberOfControlVolumes - 1, m_numberOfControlVolumes - 1]
        = 3 / m_deltax;
      mat[m_numberOfControlVolumes - 1, m_numberOfControlVolumes - 2]
        = -1 / m_deltax;
      vec[m_numberOfControlVolumes - 1, 0]
        = m_beta * m_deltax + 2 * T_R / m_deltax;

      //LUFactorization.Solve(mat, vec);
      LinearEquation.Solve(mat, vec);

      return vec;
    }

    public double posx(int i) {
      return m_deltax / 2.0 + i * m_deltax;
    }
  }
}
